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摘 要 :为 提高 流体 的 计算 效率 并 保证 结果 的 准确 性 ,利用 CUDA 编程 平台 和 GPU 强大 的 浮 点 计算 能 力 ,实现 了 基于 唱 
格 玻 尔 效 曼 方 法 的 泊 松 流 模 拟 计算 加 速 。 设 计 了 线性 寻 址 和 下 标 寻 址 2 种 不 同 寻 址 方式 ,将 这 2 种 寻 址 方式 分 别 应 用 到 
唱 格 玻 尔 效 曼 程序 的 格 点 碰撞 、 迁 徙 流动 . 安 观 量 计 算 等 步 又 中 ,并 探讨 2 种 寻 址 方式 对 程序 计算 效率 带 来 的 影响 。 同 时 
在 程序 中 使 用 统一 内 存 管 理 , 通 过 这 样 的 方式 开辟 内 存 的 变量 可 在 主机 端 和 设备 端 同时 使 用 ,简化 了 代码 复杂 度 , 同 时 降 
低 了 频繁 为 变量 开辟 内 存 带 来 的 消耗 。 使 用 Intel(R) Xeon(R) E-52620 v4 CPU,Nvidia Quadro GP100 GPU 进行 计算 ,在 

吗 性 寻 址 方法 和 下 标 寻 址 方法 中 分 别 获得 了 71 倍 和 25 倍 CPU 串 行 代码 的 加 速 比 。 
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CUDA accelerated optimization based on lattice Boltzmann method 
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sbstract: In order to improve the efficiency of fluid calculations and ensure the accuracy of the results, the CUDA program- 


‘ming platform and the powerful floating-point computing capabilities of the GPU are used to accelerate the Poisseuille flow 


_‘Sihnulation calculation based on the lattice Boltzmann method. Two different addressing methods, linear addressing and sub- 


> sEript addressing are designed, these two addressing methods are respectively applied to the lattice point collision, migration 
(Ew,and macroscopic calculation of the lattice Boltzmann program,then discuss the influence of two addressing methods on 
[ 


_He calculation efficiency of the program. At the same time, unified memory management is used in the program,and the var- 


Mables opened up in this way can be used on the host side and the device side at the same time, which simplifies the code com- 
Ovity and reduces the consumption of frequently opening up memory for variables. Using Intel(R) Xeon(R) E-52620 v4 
CPU and Nvidia Quadro GP100GPU for calculations, the linear addressing method and the subscript addressing method have 
obtained 71 times and 25 times the speedup ratio of CPU serial code respectively. 


Key words: CUDA; lattice boltzmann method; plane poiseuile flow; linear addressing; subscript addressing 


随 着 处 理 咒 时钟 频 率 极限 及 “* 功 耗 墙 ? 等 问题 的 


出 现 , 单 核 解决 方案 被 气 弃 ,基于 多 核 架 构 的 图 像 处 
理 单元 (GPU) 逐 渐 出 现在 人 们 的 视野 中 。2007 年 ， 
英 伟 达 公 司 为 GPU 增加 了 一 个 易 用 的 编程 接 
口 统一 计算 架构 (compute unified device archi- 
tecture, 简称 CUDA) 站 。 它 以 标准 C 为 基础 进行 扩 
展 , 使 程序 员 可 以 在 CC 十 十 及 Fortran 等 开发 环境 
进行 并 行程 序 编写 ,将 利用 GPU 进行 并 行 计算 的 门 


槛 大 大 降低 。 近 年 来 ,人 工 智 能 在 深度 学 习 的 不 断 推 
动 下 高 速 发 展 ,这 离 不 开 GPU 的 强大 计算 能 力 和 优 
秀 的 GPU 编程 环境 。CUDA 技术 已 经 获得 广泛 关 
注 ,并 被 应 用 到 流体 力学 .生物 计算 、 气象 分 析 、 金 融 
分 析 等 众多 领域 。 

晶 格 玻 尔 效 曼 方法 (LBM) 是 一 种 介 观 尺度 的 模 
拟 方法 ,在 流体 力学 方面 的 模拟 中 被 广泛 应 用 ,如 多 
孔 介质 流 2 .血液 流 呈 、 多 相 流 忆 淋巴 流 忆 等 。LBM 
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物理 背景 清晰 ,易于 并 行 计算 ,边界 条 件 处 理 简单 ， 
而 十 分 适用 于 大 规模 GPU 并 行 计 算 。 如 Talke 
等 中 使 用 CUDA 实现 了 LBM 多 孔 介质 流动 模型 的 
计算 加 速 。 郑 彦 奎 等 拉 实 现 了 LBM 方 腔 模型 的 计 
算 加 速 。 

鉴于 此 ,通过 线性 寻 址 和 下 标 寻 址 方法 实现 了 
LBM 模型 中 的 泊 松 流 算 例 。 将 程序 分 别 用 CPU 和 
GPU 进行 计算 ,比较 2 种 寻 址 方法 分 别 在 2 种 处 理 
单元 上 的 计算 时 间 ,并 算得 加 速 比 。 


1 GPU 与 CUDA 


基于 设计 目标 的 差异 ,GPU 在 设备 架构 上 与 
CPU 存在 很 大 区 别 。CPU 需要 很 强 的 通用 性 来 处 
理 和 从 种 不 同 的 数据 类 型 ,同时 逻辑 判断 又 会 引入 大 量 
的 窒 支 跳 转 和 中 断 处 理 。 因 此 ,CPU 中 分 布 着 多 样 
的 计算 .控制 单元 和 占据 大 量 空间 的 存储 单元 。 而 
GEC 的 设计 主要 是 用 来 解决 大 量 逻辑 上 相对 简单 的 
任 荔 ,因而 GPU 采用 了 数量 众多 的 核心 和 算术 逻辑 
单元 (ALU) ,其 数量 远 超 CPU, 通 常 CPU 的 核心 数 
不 会 超过 2 位 数 ,但 GPU 只 配备 了 简单 的 控制 逻辑 
箭 泗 化 了 的 存储 单元 ， 

〇 CUDA 编程 模型 将 CPU 作为 主机 端 ,GPU 作 
六 说 备 端 ,二 者 各 自 拥有 相互 独立 的 存储 地 址 空间 ; 
主 税 端 内 存 和 设备 端 显存 。 一 个 典型 的 CUDA 程序 
全 畦 并 行 代码 和 与 其 互补 的 串 行 代码 。 由 CPU 执 
是， 
行 蕉 码 ,并 调用 GPU 计算 密集 型 的 大 规模 数据 并 行 
计 入 代码 , 即 内 核 函 数 。 当 设备 端 开 始 执行 内 核 本 数 


时 ,设备 中 会 产生 大 量 的 线程 ,线程 是 内 核 函 数 的 基 
本 单元 , 负责 执行 内 核 函 数 的 指定 语句 站。CUDA 
程序 模型 如 图 1 所 示 ,由 一 个 内 核 启 动 所 产生 的 所 有 
线程 统称 为 一 个 线程 网 格 , 同 一 网 格 中 的 所 有 线程 共 
享 相同 的 全 局 内 存 空间 。 一 个 网 格 又 由 多 个 线程 块 
构成 ,同一 线程 块 内 的 线程 协作 可 通过 同步 和 共享 内 
存 的 方式 实现 ,不 同 块 内 的 线程 不 能 协作 。 每 个 
GPU 设备 包含 一 组 流 多 处 理 器 (SM) ,一 个 SM 又 由 
多 个 流 处 理 器 (SP) 和 一 些 其 他 硬件 组 成 。CUDA 程 
序 执行 过 程 中 会 把 线程 块 中 的 所 有 线程 以 线程 束 
(Warp, 含 32 个 线程 ) 为 执行 单位 分 配给 SM 中 的 不 
同 SP 来 进行 计算 。 


主机 端 (CPU) 设备 端 (GPU) 
线程 网 格 1 
囊 行 代码 线程 块 线程 块 线程 块 
(0,0) (1.0) (2,0) 
开行 代 三 线程 块 | 小 线程 块 | 十 线程 块 
Wo P| | 0 | ‘(1D (21) oe 
线程 网 格 2 1， 、 和 
串 行 代码 线程 | 「 线程 | 线程 块 (1.1) 
(0,0) 4 | 线程 | 线程 | 线程 | 线程 
和 民生 |, 4 (0,0) | (10) | (2.0) | (3.0) 
5 线程 线程 | | 线程 | 线程 | 线程 | 线程 
Kernel2 的 本 OD | QD | QD | GD 
中 | 线程 | 线程 | 线程 | 线程 
\| | C02) | (12) | (2.2) | (3,2) 
! 


图 1 CUDA 程序 模型 


GPU 采用 了 内 存 分 层 架 构 , 并 通过 多 级 缓存 机 
制 来 提高 读 写 效率 。CUDA 内 存 模型 如 图 2 所 示 ， 
GPU 中 主要 的 内 存 空间 包括 寄存 器 、 共 享 内 存 、 全 局 


[| 纹理 缓存 | 


友 仔 
常量 缓存 


图 2 CUDA 内 存 模型 


22060000141018 


242 桂林 电子 科技 大 学 学 报 


他 新 芭 娄 ， fbd 
Chnhlinax1Vf{ | 乍 其 


2022 年 6 月 


内 存 、 常 量 内 存 、 纹 理 内 存 和 本 地 内 存 。 每 个 SM 都 
拥有 独立 的 一 级 缓存 ,所 有 SM 共享 二 级 缓存 。 
GPU 致力 于 为 每 个 线程 分 配 真实 的 寄存 融 , 当 寄存 
器 使 用 到 达 上 限时 ,编译 器 会 将 数据 放 在 片 外 的 本 地 
内 存 中 。 共 享 内 存 是 可 受 程 序 员 控制 的 一 级 缓存 ,每 
es 中 的 一 级 缓存 与 共享 内 存 公用 同一 个 内 存 
股 ” , 它 仅 对 正在 执行 的 线程 块 中 的 每 个 线程 可 见 。 


2 晶 格 玻 尔 兹 曼 方 法 (LBMD) 


McNamara 等 :用 单 粒 子 分 布 函数 f; 取代 了 
格子 气 细 胞 自动 机 中 的 布尔 变量 ,其 LBM 演化 方 
程 为 

fednt to) — frst)=0(f)s (1) 
其 中 .i 为 微观 速度 的 索引 ,对 于 D2Q9 模型 ,i 取 值 
EY9; fi(x ,it) 为 在 x 位 置 1 时 刻 具 有 e; 速度 的 粒子 
的 分 布 函 数 ; 5, 为 时 间 增 量 ; Q(f;) 为 碰撞 因子 , 表 
下 三 撞 对 于 /， 的 影响 。 文 献 L11-13] 采 用 单 凶 豫 时 
问 邹 痊 代 碰撞 因子 项 ， 继而 ,LBM 方程 可 写 为 


1 _. 
ot +60— fix = ff) ， 


/@Q+e8 


(2) 


其 昌 . /, 为 局 域 平衡 分 布 函 数 ; r 为 弛 随时 间 。 求 
时 了 分布 数 后 , 则 宏观 密度 ,流体 的 流速 分 别 为 
0=2 fi=2 1", C3) 
i Deif 2 
Ee u 一 一 (4) 
一 O O 


玫 则 因 3 所 示 的 D2Q9 模型 (D 指 维 度 ,Q 指 粒子 运 
动 月 向 总 数 ) 模 型 " 呈 进 行 计算 ,其 平衡 分 布 函数 具体 
定义 为 


了 一 OO | | (e€ju) | 
C 


(Cej&D) 
2c 


| | 
ee 上 te 
7， 4 ‘Cg 


图 3 D2Q9 模型 


其 中 :c 为 基准 速度 ,通常 情况 下 取 1; w; 为 各 方向 的 
权重 系数 。 当 i 二 0 时 , w,; 二 4/9; 当 1 二 1,2,3,4 时 ， 


w; 二 1/9; 当 i 二 5,6,7,8 时 , w; 二 1/36。e@; 的 形式 为 

= 0 1 一 1 —1 "| 
my 0 =L1 1T TT = 
3 算法 实现 及 优化 


一 个 典型 的 LBM 计算 案例 主要 分 为 3 个 步 又: 
了 ) 声 明 变 量 并 进行 分 布 函 数 的 初始 化 ;2) 进 行 格 点 的 
碰撞 、 迁 徙 流动 和 边界 处 理 , 并 进行 宏观 量 的 计算 ; 
3) 进 行 稳定 性 条 件 判 断 ,决定 是 否 结束 步骤 2) 的 迭 
代 。 由 文献 [15j 知 ,步骤 2) 的 迭代 计算 占 整 个 计算 
时 间 的 98%% ,因此 该 步骤 应 在 GPU 上 进行 并 行 加 
速 。 本 计算 以 淋巴 管 为 物理 背景 ,将 淋巴 管内 流动 的 
淋巴 液 视 作 等 截面 圆 形 管道 内 的 泊 松 流 ,进行 模拟 计 
算 。 尝 试用 线性 寻 址 和 下 标 寻 址 2 种 不 同 的 寻 址 方 
法 进行 计算 , 比较 这 2 种 寻 址 方法 分 别 在 GPU 与 
CPU 上 执行 计算 时 间 的 差异 。 


3.1 线性 寻 址 程序 设计 


LBM 的 基本 变量 是 分 布 函数 ,进而 由 分 布 函 数 
计算 求 得 格子 点 的 密度 和 zy 方向 上 的 速度 等 宏观 
量 。 在 串 行程 序 的 D2Q9 模型 中 ,分 布 函 数 通 常 被 设 
为 三 维 数组 dfLijUjjLkgj 并 存储 在 主 存 中 ,表示 整 
个 线程 网 格 中 格 点 的 坐标 (0 二 i 二 width, 0 入 7 三 
height; width、height 均 为 固定 值 ,分 别 表示 线程 网 格 
的 宽度 和 高 度 ) ,k 为 格 点 运动 方向 数 ,上 一 0,1，…,8。 

在 设备 端 声明 指针 double x df ,调用 cudaMal- 
locManaged( ) 函数 在 设备 端 开辟 大 小 为 width x 
height x 9 x sizeof (double) 的 一 维 线性 内 存 空间 ,用 
于 存放 每 个 格 点 上 9 个 方向 的 分 布 函数 ,并 将 在 显存 
上 获得 的 内 存 空 间 首 地 址 赋值 给 df。 使 用 cudaMal- 
locManaged() 函 数 开 辟 的 存储 空间 ,无 论 是 在 串 行 
代码 中 还 是 并 行 代码 中 ,都 可 使 用 这 块 内 存 ,因此 只 

需 定义 一 个 指针 即 可 在 主机 和 设备 端 通用 。 继 续 开 
辟 多 个 设备 端 内 存 空间 ,用 于 存放 其 他 计算 变量 的 数 
据 。 指 针 double * dd 指向 存储 格子 点 密度 的 内 存 
空间 ,指针 double * dvx、dvy 分 别 指向 存放 之 和 xy 
方向 上 速度 的 内 存 空 间 。 

基于 上 述 在 全 局 内 存 中 使 用 三 维 数组 来 储存 分 
布 函 数 及 其 他 变量 的 方式 ,设计 了 3 个 在 CPU 端 执 
行 的 collide() stream() 、calculate() 也 数 ,分 别 代 表 
格 点 的 碰撞 、 迁 徙 流动 和 宏观 量 的 计算 等 步骤 ,并 将 
其 命名 为 方案 一 。 方 案 二 则 是 将 迭代 计算 过 程 放 在 
GPU 端 并 行 执行 ,对 应 设计 了 addKernelCollide()、 
addKernelCopy( )、addKernelStream ( )、 addKernel- 
Calculate() 四 个 内 核 函 数 来 实现 。 相 较 于 方案 一 中 
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执行 各 个 步骤 都 需要 骨 套 for 循环 来 遍历 读 写 每 个 
格 点 上 的 数据 ,方案 二 则 通过 内 核 浮 数 用 GPU 映射 
出 大 量 线程 ,同时 对 每 个 格 点 数据 进行 操作 。 

首先 ,调用 cudaMallocManaged() 函数 ,在 GPU 
端 为 分 布 函数 、 密 度 等 变量 开辟 内 存 空间 ,再 定义 一 
个 parameter 结构 体 , 为 之 后 格子 点 迁徙 流动 做 铺 
垫 。 对 4 个 内 核 函 数 进 行 如 下 讨论 : 

1) 碰 撞 步 ;启动 addKernelCollide 志 < 二 到 grid， 
block>>>>>(df,iSol,…,dev_p) 内 核 函 数 ,其 中 grid 
为 线程 网 格 中 线程 块 数 ,将 其 设 为 width; block 为 每 
个 线程 块 中 的 线程 数 , 设 为 9 * height;iSol 为 用 于 判 
上 断 格 点 是 否 位 于 边界 上 的 变量 ;dev_p 为 parameter 
结构 体 变量 。 声 明 变 量 i、j、k 并 赋值 ,i=blockIdx. 
x= threadIdx. y， k 二 threadIdx. x。 在 函数 体内 分 
别 渤 算出 线程 访问 分 布 函数 一 维 数组 的 下 标 索 引 
idssadr(i,j,k) 一 k 十 (十 ix height) x* 9 和 其 他 变量 
的 下 标 索引 ind 二 adr(i,j) 二 十 1x 1x height。 通 过 
一 次 站 判断 ,保留 格子 边界 以 外 的 所 有 格子 点 进行 
计 血 ,每 个 线程 利用 式 (2) 计 算 对 应 格子 点 的 分 布 
函数 . 

2 df[id]=df[id]— (df[idj— 
四 feqk,dd[ind],dvx[ind],dvy[ind]))/Tau, 
f 人 0 四 为 平衡 分 布 函 数 ,Tau 为 弛 驳 时 间 。 

> 2) 流 动迁 徙 步 : 启 动 addKernelCopy 二 二 过 grid， 
ble 和 二 二 (df,dfbak) 内 核 函数 ,声明 变量 并 赋值 , 计 
算 饪 线程 访问 的 下 标 索 引 , 方 法 同上 。 将 df[idj 中 的 
数据 传递 给 dfbak| id ]|。 启 动 addKernelStream 二 二 
-<< 和 tid ,block 这 之 (df,dfbak,iSol,…) 内 核 函 数 ,将 
dfbak[idj 的 值 赋 给 df[idj, 此 时 dfbak[id] 中 

id 一 kk 十 LG 一 dev_PrjyLk]) 十 (Ci 一 
dev_Prjx[ k |) * height) | * 9。 

3) 计 算 宏 观 量 步 :启动 addKernelCalculate 二 二 
二 grid, block 放 放 之 (df,iSol,…) 内 核 函数 ,声明 变量 
并 放 入 共享 内 存 中 。 

__shared ”double sdf[9|];//9 个 分 布 函 数 
_ shared_ ”double ddd; // 格 子 点 密度 
__shared double vx;//Zx 方向 速度 
__shared double vy;//y 方向 速度 
for (k=0;k<=9;k 十 十 ) 

( 

sdf[k |=dfladr(i,j, k) |; 

} 
利用 sdf[ kj 计算 得 到 ddd、vx、vy 的 值 ,最 后 将 共享 
内 存 中 的 数据 传 回 主 存 中 。 


3. 2 ”下 标 寻 址 程序 设计 


基于 之 前 的 线性 寻 址 方法 对 程序 做 以 下 修改 。 

1) 声 明 指针 double x df0、double xx dfl ,double 
xx xx df, 在 GPU 上 开辟 3 个 内 存 空 间 ,并 将 首 地 址 分 
别 赋值 给 dfo .dfl df, 具体 操作 如 下 : 
cudaMallocManaged(void(x* x* ) &.df0, width x height 
x 9 x sizeof (double)); 
cudaMallocManaged(void(x * ) &.df1 ,width x height 
x* sizeof (double * ) ); 
for (i 二 0;1 达 width;i 十 十 ) 

( 
for(j 二 0;] 过 height;j 十 十 ) 
{dfl[ix height+j|= &df0[ (i x height 十 j) 
x* 9];) 
} 
cudaMallocManaged( void (x * ) &df, width * height 
x sizeof (double * * )); 
for (i 二 0;1 过 width;i 十 十 ) 
{df[i]= &.dfl[ix height |;} 

2) 声 明 指 针 double x* dd0、double *x dd, 在 
GPU 上 开辟 2 次 内 存 空 间 , 并 将 首 地 址 分 别 赋 值 给 
dd0、dd, 具 体操 作 如 下 : 
cudaMallocManaged (void (xx ) &dd0, width * 
height x Sizeof(double) ) ; 
cudaMallocManaged (void (x * ) &.dd, width * sizeof 
(double * )); 
for (i 二 0;1 达 width;i 十 十 ) 

{dd[i]= &.ddO[ix height |;} 

3) 其 他 宏观 量 修 改 方法 同上 。 

在 各 计算 步骤 中 ,用 多 维 数组 表示 不 同 格 点 位 置 
上 的 分 布 函数 和 其 他 宏观 量 , 形 如 df[ij[jj[Lk]、dd[i] 
jdvx[ ij[j]\dvy[LijLjj。 将 方案 一 、 二 按 上 述 内 容 
进行 修改 ,并 分 别 命 名 为 方案 三 、 四 。 


4 计算 结果 分 析 


采用 含有 8 个 Nvidia Quadro GP100 显卡 集群 
的 服务 器 完成 计算 ,GP100 显卡 拥有 3 584 个 CUDA 
并 行 计算 处 理 核 心 ,处 理 双 精 度 浮 点 数 的 能 力 为 5. 2 
TFlop/s, 搭 配 16 GiB HBM2 显存 ,理论 带宽 高 达 
717 GiB/s。 同 时 服务 器 搭载 了 Intel(R) Xeon(R) 
CPU E-52620 v4 处 理 器 ,核心 数 为 8 个 , 主 频 为 2.1 
GHz, 编译 环境 为 CUDA10.0, 运行 环境 为 Linux 
Ubuntu。 

以 平面 泊 松 流 为 算 例 , 在 保证 计算 结果 正确 性 的 
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